function delta = computeGradient(i,robot,r,dr_dh,points_xy,mesh_size,R,d_center)
    s_i = robot(i,:); state_dim = size(robot,2); r_num = size(robot,1);
    % search for the neighboring robots
%     neighborRobots = vecnorm(s_i-robot,2,2) < r(i)+r;
%     neighborRobots(i) = 0;
%     neighborRobots = find(neighborRobots==1);
%     delta = zeros(3,1);
    
    % refine the integral area to the coverage of s_i
    integral_index = ( vecnorm([points_xy zeros(size(points_xy,1),1)] - s_i,2,2)<=r(i) );
    integral_xy = points_xy( integral_index,: );
    integral_dim = size(integral_xy,1);
    % repeat the process in computeObj_vecform to obtain the PDOP
    % the derivative in coverage area
    delta_area = computeGradientArea(i,robot,integral_xy,integral_index,integral_dim,r,R,mesh_size);
    
    % the derivative in boundary
    
    theta = (0:pi/100:2*pi)';
    theta_size = pi/100;
   
    integral_xy = [r(i)*cos(theta)+s_i(1) r(i)*sin(theta)+s_i(2)];
    integral_dim = size(integral_xy,1);
    delta_boundary = computeGradientBoundary(i,robot,theta,integral_xy,integral_dim,r,theta_size,dr_dh,d_center);

    delta = delta_area + delta_boundary;
end

function delta_boundary = computeGradientBoundary(r_index,robot,theta,integral_xy,integral_dim,r,theta_size,dr_dh,d_center)
[r_num,state_dim] = size(robot);
delta_boundary = zeros(state_dim,1);    
count_num = zeros(integral_dim,r_num);
    for i = 1:r_num
        if i == r_index
            count_num(:,i) = ones(integral_dim,1);
            continue;
        end
        posi = integral_xy - robot(i,1:2);
        dis = vecnorm(posi,2,2);
        count_num(:,i) = (dis<=r(i));    
    end
    zx = zeros(integral_dim,1); zy = zx; zz = zx;
    
    r1 = sum(count_num,2) == 1;    
    k1 = find(r1==1);
    for j = 1:size(k1)
        kk = k1(j);
        pdop = 0.05;
        zx(kk) = pdop*r(r_index)*cos( theta(kk) );
        zy(kk) = pdop*r(r_index)*sin( theta(kk) );
        zz(kk) = pdop*r(r_index)*dr_dh(r_index);
    end

    r1 = sum(count_num,2) == 2; 
    k1 = find(r1==1);
    for j = 1:size(k1)
        kk = k1(j);
        pdop = 0.1;
        zx(kk) = pdop*r(r_index)*cos( theta(kk) );
        zy(kk) = pdop*r(r_index)*sin( theta(kk) );
        zz(kk) = pdop*r(r_index)*dr_dh(r_index);
    end
    
    r1 = sum(count_num,2) >=3; r3 = kron(r1,ones(r_num,1)); 
    r2 = reshape(count_num',numel(count_num),1);  r2 = r2&r3;
    
    robot_rel = kron(ones(integral_dim,1),robot);
    points = [kron(integral_xy,ones(r_num,1)) zeros(r_num*integral_dim,1)];

    G = (points - robot_rel);
    G = G.*(r2)./vecnorm(G,2,2);

    k1 = find(r1==1);
    for j = 1:size(k1)
        kk = k1(j);
        Gi = G(r_num*(kk-1)+1:r_num*(kk-1)+r_num,:);
        p = trace(inv(Gi'*Gi)); p = sigmoid(p);
        zx(kk) = p*r(r_index)*cos( theta(kk) );
        zy(kk) = p*r(r_index)*sin( theta(kk) );
        zz(kk) = p*r(r_index)*dr_dh(r_index);
    end
    center = d_center;
    R_boundary = densityFunc(center,integral_xy);
    delta_boundary(1) = sum(R_boundary.*zx*theta_size,'all');
    delta_boundary(2) = sum(R_boundary.*zy*theta_size,'all');
    delta_boundary(3) = sum(R_boundary.*zz*theta_size,'all');
end

function delta_area = computeGradientArea(r_index,robot,integral_xy,integral_index,integral_dim,r,R,mesh_size)
    delta_area = zeros(size(robot,2),1); s_i = robot(r_index,:);
    zx = zeros(integral_dim,1);zy = zx; zz = zx;
    if size(robot,2) == 1
        robot = reshape(robot,numel(robot)/3,3);
    end
    [r_num,state_dim] = size(robot);  
    count_num = zeros(integral_dim,r_num);
    for i = 1:r_num
        if i == r_index
            count_num(:,i) = ones(integral_dim,1);
            continue;
        end
        posi = integral_xy - robot(i,1:2);
        dis = vecnorm(posi,2,2);
        count_num(:,i) = (dis< r(i));    
    end
   
    r1 = sum(count_num,2) >=3; r3 = kron(r1,ones(r_num,1)); 
    r2 = reshape(count_num',numel(count_num),1);  r2 = r2&r3;
    
    robot_rel = kron(ones(integral_dim,1),robot);
    points = [kron(integral_xy,ones(r_num,1)) zeros(r_num*integral_dim,1)];

    G = (points - robot_rel);
    G = G.*(r2)./vecnorm(G,2,2);
    
    k1 = find(r1==1);
    for j = 1:size(k1)
        k = k1(j);
        q = integral_xy(k,:);
        Gi = G(r_num*(k-1)+1:r_num*(k-1)+r_num,:);
        [~,ds_dpdop] = sigmoid( pdop(Gi) );
        [Mx,My,Mz] = d_gTg_d_s(q,s_i);        
        zx(k) = -ds_dpdop*trace(inv(Gi'*Gi)*Mx*inv(Gi'*Gi));     % d pdop/ d s_i
        zy(k) = -ds_dpdop*trace(inv(Gi'*Gi)*My*inv(Gi'*Gi));
        zz(k) = -ds_dpdop*trace(inv(Gi'*Gi)*Mz*inv(Gi'*Gi));
    end
    R = reshape(R,numel(R),1);
    integral_R = R(integral_index);
    delta_area(1) = sum(integral_R.*zx*mesh_size^2,'all');
    delta_area(2) = sum(integral_R.*zy*mesh_size^2,'all');
    delta_area(3) = sum(integral_R.*zz*mesh_size^2,'all');
end

function [Mx,My,Mz] = d_gTg_d_s(q,s_i)
% compute the derivative of d (gi'gt)/ (d s_x)
    d = norm( [q 0] - s_i,2);
    Mx = (1/d^2) * [ 
        2*( s_i(1) -q(1) )      ( s_i(2) - q(2) )    s_i(3);
        ( s_i(2) - q(2))                    0                  0;
        s_i(3)                                 0                  0;
    ];
    My = (1/d^2)* [ 
        0                               (s_i(1) - q(1))             0;
        (s_i(1) - q(1))           2*(s_i(2)-q(2))           s_i(3);
        0                                   s_i(3)                      0;
    ];
    Mz = (1/d^2)* [
        0                                   0                      (s_i(1) - q(1));
        0                                   0                      (s_i(2) - q(2));
        (s_i(1) - q(1))        (s_i(2) - q(2))                  2*s_i(3);
    ];
end

